<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.1 plus MathML 2.0//EN" "http://www.w3.org/Math/DTD/mathml2/xhtml-math11-f.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" lang="en-US">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8"/>
<meta name="viewport" content="width=device-width, initial-scale=0.4"/>
<meta name="google" content="notranslate" />
<link rel="canonical" href="https://sleef.org/additional.xhtml" />
<link rel="icon" href="favicon.png" />
<link rel="stylesheet" type="text/css" href="texlike.css"/>
<link rel="stylesheet" type="text/css" href="//fonts.googleapis.com/css?family=Ubuntu" />
<link rel="stylesheet" type="text/css" href="sleef.css"/>
<title>SLEEF - Additional Notes</title>
</head>
<body translate="no" class="notranslate">
<h1>SLEEF Documentation - Additional Notes</h1>

<h2>Table of contents</h2>

<ul class="none" style="font-family: arial, sansserif; padding-left: 0.5cm;">
  <li><a class="underlined" href="index.xhtml">Introduction</a></li>
  <li><a class="underlined" href="compile.xhtml">Compiling and installing the library</a></li>
  <li><a class="underlined" href="purec.xhtml">Math library reference</a></li>
  <li><a class="underlined" href="dft.xhtml">DFT library reference</a></li>
  <li><a class="underlined" href="misc.xhtml">Other tools included in the package</a></li>
  <li><a class="underlined" href="benchmark.xhtml">Benchmark results</a></li>
  <li>&nbsp;</li>
  <li><a class="underlined" href="additional.xhtml">Additional notes</a></li>
    <ul class="disc">
      <li><a href="additional.xhtml#faq">Frequently asked questions</a></li>
      <li><a href="additional.xhtml#gnuabi">About the GNUABI version of the library</a></li>
      <li><a href="additional.xhtml#dispatcher">How the dispatcher works</a></li>
      <li><a href="additional.xhtml#ulp">ULP, gradual underflow and flush-to-zero mode</a></li>
      <li><a href="additional.xhtml#paynehanek">Explanatory source code for the modified Payne Hanek reduction method</a></li>
      <li><a href="additional.xhtml#logo">About the logo</a></li>
    </ul>
</ul>

<h2 id="faq">Frequently asked questions</h2>

<p class="noindent">
  <b>Q1:</b> Is the scalar functions in SLEEF faster than the
  corresponding functions in the standard C library?
</p>

<br/>

<p class="noindent">
  <b>A1:</b> No. Todays standard C libraries are very well optimized,
  and there is small room for further optimization. The reason why
  SLEEF is fast is that it carries out computation directly on SIMD
  registers and ALUs. This is not simple as it sounds, because
  conditional branches have to be eliminated in order to take full
  advantage of SIMD computation. If the algorithm requires conditional
  branches according to the argument, it must prepare for the case
  where the elements in the input vector contain both values that
  would make a branch happen and not happen. This would spoil the
  advantage of SIMD computation, because each element in a vector
  would require a different code path.
</p>

<br/>
<br/>

<p class="noindent">
  <b>Q2:</b> Do the trigonometric functions (e.g. sin) in SLEEF return
  correct values for the whole range of inputs?
</p>

<br/>

<p class="noindent">
  <b>A2:</b> Yes. SLEEF does implement a vectorized version of Payne
  Hanek range reduction, and all the trigonometric functions return
  a correct value with the specified accuracy.
</p>

<h2 id="gnuabi">About the GNUABI version of the library</h2>

<p class="noindent">
  The GNUABI version of the library (libsleefgnuabi.so) is built for
  x86 and aarch64 architectectures. This library provides an API
  compatible with <a class="underlined"
  href="https://sourceware.org/glibc/wiki/libmvec">libmvec</a> in
  glibc, and the API comforms to the <a class="underlined"
  href="https://sourceware.org/glibc/wiki/libmvec?action=AttachFile&amp;do=view&amp;target=VectorABI.txt">x86
  vector ABI</a>, <a class="underlined"
  href="https://developer.arm.com/docs/101129/latest">AArch64 vector
  ABI</a> and <a class="underlined"
  href="https://github.com/power8-abi-doc/vector-function-abi/">Power
  Vector ABI</a>. This library is built and installed by default, and
  certain compilers call the functions in this library.
</p>


<h2 id="dispatcher">How the dispatchers work</h2>

<p class="noindent">
  The dispatchers in SLEEF are designed to have very low overhead. This
  overhead is so small and cannot be observed by microbenchmarking.
</p>

<p>
  Fig. 7.1 shows a simplified code of our dispatcher. There is only
  one exported function <b class="func">mainFunc</b>. When
  <b class="func">mainFunc</b> is called for the first
  time, <b class="func">dispatcherMain</b> is called internally,
  since <i class="var">funcPtr</i> is initialized to the pointer to
  <b class="func">dispatcherMain</b> (line 14). It then detects if the
  CPU supports SSE 4.1 (line 7), and
  rewrites <i class="var">funcPtr</i> to a pointer to the function
  that utilizes SSE 4.1 or SSE 2, depending on the result of CPU
  feature detection (line 10).  When
  <b class="func">mainFunc</b> is called for the second time, it does
  not execute the
  <b class="func">dispatcherMain</b>. It just executes the function
  pointed by the pointer stored in <i class="var">funcPtr</i> during
  the execution of
  <b class="func">dispatcherMain</b>.
</p>

<p>
  There are advantages in our dispatcher. The first advantage is that
  it does not require any compiler-specific extension. The second
  advantage is simplicity. There are only 18 lines of simple
  code. Since the dispatchers are completely separated for each
  function, there is not much room for bugs to get in.
</p>

<p>
  The third advantage is low overhead. You might think that the
  overhead is one function call including execution of the prologue
  and the epilogue. However, modern compilers are smart enough to
  eliminate redundant execution of the prologue, epilogue and return
  instruction. The actual overhead is just one jmp instruction, which
  has very small overhead since it is not conditional. This overhead
  is likely hidden by out-of-order execution.
</p>

<p>
  The fourth advantage is thread safety. There is only one variable
  shared among threads, which is <i class="var">funcPtr</i>. There are
  only two possible values for this pointer variable. The first value
  is the pointer to the <b class="func">dispatcherMain</b>, and the
  second value is the pointer to either <b class="func">funcSSE2</b>
  or <b class="func">funcSSE4</b>, depending on the availability of
  extensions. Once <i class="var">funcPtr</i> is substituted with the
  pointer to <b class="func">funcSSE2</b>
  or <b class="func">funcSSE4</b>, it will not be changed in the
  future. It should be easy to confirm that the code works in all the
  cases.
</p>


<pre class="code">
<code>static double (*funcPtr)(double arg);</code>
<code></code>
<code>static double dispatcherMain(double arg) {</code>
<code>    double (*p)(double arg) = funcSSE2;</code>
<code></code>
<code>#if the compiler supports SSE4.1</code>
<code>    if (SSE4.1 is available on the CPU) p = funcSSE4;</code>
<code>#endif</code>
<code></code>
<code>    funcPtr = p;</code>
<code>    return (*funcPtr)(arg);</code>
<code>}</code>
<code></code>
<code>static double (*funcPtr)(double arg) = dispatcherMain;</code>
<code></code>
<code>double mainFunc(double arg) {</code>
<code>    return (*funcPtr)(arg);</code>
<code>}</code>
</pre>
<p style="text-align:center; margin-bottom: 1.0cm;">
  Fig. 7.1: Simplified code of our dispatcher
</p>


<h2 id="ulp">ULP, gradual underflow and flush-to-zero mode</h2>

<p class="noindent">
  ULP stands for "unit in the last place", which is sometimes used for
  representing accuracy of calculation. 1 ULP is the distance between
  the two closest floating point number, which depends on the exponent
  of the FP number. The accuracy of calculation by reputable math
  libraries is usually between 0.5 and 1 ULP. Here, the accuracy means
  the largest error of calculation. SLEEF math library provides
  multiple accuracy choices for most of the math functions. Many
  functions have 3.5-ULP and 1-ULP versions, and 3.5-ULP versions are
  faster than 1-ULP versions. If you care more about execution speed
  than accuracy, it is advised to use the 3.5-ULP versions along with
  -ffast-math or "unsafe math optimization" options for the compiler.
</p>

<p>
  Note that 3.5 ULPs of error is small enough in many applications. If
  you do not manage the error of computation by carefully ordering
  floating point operations in your code, you would easily have that
  amount of error in the computation results.
</p>

<p>
  In IEEE 754 standard, underflow does not happen abruptly when the
  exponent becomes zero. Instead, when a number to be represented is
  smaller than a certain value, a denormal number is produced which
  has less precision. This is sometimes called gradual underflow. On
  some processor implementation, a flush-to-zero mode is used since it
  is easier to implement by hardware. In flush-to-zero mode, numbers
  smaller than the smallest normalized number are replaced with
  zero. FP operations are not IEEE-754 conformant if a flush-to-zero
  mode is used. A flush-to-zero mode influences the accuracy of
  calculation in some cases. The smallest normalized precision number
  can be referred with DBL_MIN for double precision, and FLT_MIN for
  single precision. The naming of these macros is a little bit
  confusing because DBL_MIN is not the smallest double precision
  number.
</p>

<p>
  You can see known maximum errors in math functions in glibc
  at <a class="underlined"
  href="http://www.gnu.org/software/libc/manual/html_node/Errors-in-Math-Functions.html">
  this page.</a>
</p>


<h2 id="paynehanek">Explanatory source code for our modified Payne Hanek reduction method</h2>

<p class="noindent">
  In order to evaluate a trigonometric function with a large argument,
  an argument reduction method is used to find an FP remainder of
  dividing the argument <i class="var">x</i> by &pi;. We devised a
  variation of the Payne-Hanek argument reduction method which is
  suitable for vector computation. Fig. 7.2
  shows <a class="underlined" href="ph.c">an explanatory source
  code</a> for this method. See <a class="underlined"
  href="http://dx.doi.org/10.1109/TPDS.2019.2960333">our paper</a> for
  the details.
</p>

<pre class="code">
<code>#include &#60;stdio.h&#62;</code>
<code>#include &#60;stdlib.h&#62;</code>
<code>#include &#60;math.h&#62;</code>
<code>#include &#60;mpfr.h&#62;</code>
<code></code>
<code>typedef struct { double x, y; } double2;</code>
<code>double2 dd(double d) { double2 r = { d, 0 }; return r; }</code>
<code>int64_t d2i(double d) { union { double f; int64_t i; } tmp = {.f = d }; return tmp.i; }</code>
<code>double i2d(int64_t i) { union { double f; int64_t i; } tmp = {.i = i }; return tmp.f; }</code>
<code>double upper(double d) { return i2d(d2i(d) &#38; 0xfffffffff8000000LL); }</code>
<code>double clearlsb(double d) { return i2d(d2i(d) &#38; 0xfffffffffffffffeLL); }</code>
<code></code>
<code>double2 ddrenormalize(double2 t) {</code>
<code>  double2 s = dd(t.x + t.y);</code>
<code>  s.y = t.x - s.x + t.y;</code>
<code>  return s;</code>
<code>}</code>
<code></code>
<code>double2 ddadd(double2 x, double2 y) {</code>
<code>  double2 r = dd(x.x + y.x);</code>
<code>  double v = r.x - x.x;</code>
<code>  r.y = (x.x - (r.x - v)) + (y.x - v) + (x.y + y.y);</code>
<code>  return r;</code>
<code>}</code>
<code></code>
<code>double2 ddmul(double x, double y) {</code>
<code>  double2 r = dd(x * y);</code>
<code>  r.y = fma(x, y, -r.x);</code>
<code>  return r;</code>
<code>}</code>
<code></code>
<code>double2 ddmul2(double2 x, double2 y) {</code>
<code>  double2 r = ddmul(x.x, y.x);</code>
<code>  r.y += x.x * y.y + x.y * y.x;</code>
<code>  return r;</code>
<code>}</code>
<code></code>
<code>// This function computes remainder(a, PI/2)</code>
<code>double2 modifiedPayneHanek(double a) {</code>
<code>  double table[4];</code>
<code>  int scale = fabs(a) > 1e+200 ? -128 : 0;</code>
<code>  a = ldexp(a, scale);</code>
<code></code>
<code>  // Table genration</code>
<code></code>
<code>  mpfr_set_default_prec(2048);</code>
<code>  mpfr_t pi, m;</code>
<code>  mpfr_inits(pi, m, NULL);</code>
<code>  mpfr_const_pi(pi, GMP_RNDN);</code>
<code></code>
<code>  mpfr_d_div(m, 2, pi, GMP_RNDN);</code>
<code>  mpfr_set_exp(m, mpfr_get_exp(m) + (ilogb(a) - 53 - scale));</code>
<code>  mpfr_frac(m, m, GMP_RNDN);</code>
<code>  mpfr_set_exp(m, mpfr_get_exp(m) - (ilogb(a) - 53));</code>
<code></code>
<code>  for(int i=0;i&#60;4;i++) {</code>
<code>    table[i] = clearlsb(mpfr_get_d(m, GMP_RNDN));</code>
<code>    mpfr_sub_d(m, m, table[i], GMP_RNDN);</code>
<code>  }</code>
<code></code>
<code>  mpfr_clears(pi, m, NULL);</code>
<code></code>
<code>  // Main computation</code>
<code></code>
<code>  double2 x = dd(0);</code>
<code>  for(int i=0;i&#60;4;i++) {</code>
<code>    x = ddadd(x, ddmul(a, table[i]));</code>
<code>    x.x = x.x - round(x.x);</code>
<code>    x = ddrenormalize(x);</code>
<code>  }</code>
<code></code>
<code>  double2 pio2 = { 3.141592653589793*0.5, 1.2246467991473532e-16*0.5 };</code>
<code>  x = ddmul2(x, pio2);</code>
<code>  return fabs(a) &#60; 0.785398163397448279 ? dd(a) : x;</code>
<code>}</code>
</pre>
<p style="text-align:center; margin-bottom: 1.0cm;">
  <a href="ph.c">Fig. 7.2: Explanatory source code for our modified Payne Hanek reduction method</a>
</p>



<h2 id="logo">About the logo</h2>

<p>
  It is a soup ladle. A sleef means a soup ladle in Dutch.
</p>

<br/>

<p style="text-align:center; margin-top:1cm;">
  <a class="nothing" href="sleeflogo3.svg">
    <img src="sleeflogo3.png" alt="logo" width="40%" height="40%" />
  </a>
  <br />
  Fig. 7.2: SLEEF logo
</p>

<p class="footer">
  Copyright &copy; <!--YEAR--> SLEEF Project.<br/>
  SLEEF is open-source software and is distributed under the Boost Software License, Version 1.0.
</p>

<!--TEST-->

</body>
</html>
